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We use a simple and efficient computer model to investigate the physical properties of bilayer 
membranes. The amphiphilic molecules are modeled as short rigid trimers with finite range pair 
interactions between them. The pair potentials have been designed to mimic the hydrophobic 
interactions, and to allow the simulation of the membranes without the embedding solvent as if 
the membrane is in vacuum. We find that upon decreasing the area density of the molecules the 
membrane undergoes a solid-fluid phase transition, where in the fluid phase the molecules can diffuse 
within the membrane plane. The surface tension and the bending modulus of the fluid membranes 
are extracted from the analysis of the spectrum of thermal undulations. At low area densities we 
observe the formation of pores in the membrane through which molecules can diffuse from one layer 
to the other. The appearance of the pores is explained using a simple model relating it to the area 
dependence of the free energy. 

PACS numbers: 



I. INTRODUCTION 

When amphiphilic molecules such as lipids are brought 
into contact with water they tend to arrange so as to 
shield their "oily" hydrocarbon tail from the aqueous en- 
vironment while exposing their hydrophilic head to the 
water. One of the simplest structures formed in this way 
is of a bilayer membrane - a double sheet of surfactants 
separating two aqueous phases m. Bilayer membranes 
are common in biological systems ■ Living cells are sep- 
arated from their extra-cellular surroundings by plasma 
membranes that control the transport of material into 
and out of the cell 0, 0| • Most biological membranes are 
found in the fluid phase where the lipids comprising the 
bilayer can diffuse freely in the membrane plane. Another 
characteristic feature of lipid bilayers is their high flexibil- 
ity which allows for large thermally-excited undulations 
pl, |6j . The fluidity and low rigidity of membranes are im- 
portant for many of their biological properties, such as 
their ability to change their shape easily and the possi- 
bility of proteins to insert themselves into the membrane 

The thickness of membranes is comparable to the size 
of the constituting surfactant molecules (typically on 
the nanometer scale), while their lateral extension can 
greatly exceed their thickness and reach up to several mi- 
crometers. Consequently, coarse-grain phenomenological 
models, such as Ginzburg-Landau free energy function- 
als [1| or the effective surface Hamiltonian IE IS 13 ' 
have been used in order to study the physical properties 
of membranes, as well as of other interfaces (like surfac- 
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tants monolayer in microemulsions or vapor-fluid inter- 
faces). In those theories the bilayer membrane is treated 
as a smooth continuous surface, and its elastic energy is 
related to the membrane area and the local curvatures. 
These theories have been very successful in descri bing th e 
shape and phase diagrams of bilayer membranes llll | . 

Phenomenological models describe the mesoscopic 
physical behavior of interfaces and membranes, but do 
not allow one to approach these systems on the molecular 
level. Many theories have been developed in an attempt 
to understand how the mesoscopic behavior emerges from 
the microscopic entities and the interactions between 
them. These theories include lattice "Ising-like" models 
[a, molecular theories of the hydrocarbon chain packing 
theories including the effect of electrostatic inter- 
actions and density functional theories 0|. The 
most microscopic detailed approach is employed in some 
computer simulations where the amphiphiles and water, 
and the interactions between them are modeled explic- 
itly in full detail [3. Since these simulations require 
an enormously large computing time, they are restricted 
to fairly small systems consisting of 50-200 amphiphiles, 
and can be utilized to investigate phenomena occurring 
on short time scales of a few nanoseconds. In order to 
study mesoscale phenomena it is therefore necessary to 
dispense with some of the microscopic details in the simu- 
lations and use simplified models [13 • A number of such 
simplified computer models have been devised by several 
groups. In these models the structure of the surfactant 
molecules is represented in a "coarse-grained" manner 
where a number of atoms are grouped together into a sin- 
gle site. The first level of coarse-graining is obtained by 
replacing the water molecules and the CH2 groups of the 
hydrocarbon chain by unified atoms |i3illall3- This can 
reduce the number of atoms per lipids to about 50. Much 
more simplified models, in which the amphiphiles consists 
of only 5-10 atoms, were also presented [23,|^|22|- In 
these latter models the electrostatic potentials are usu- 
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ally ignored and the potentials of the chemical bonds are 
greatly simplified. At this level of simplification it is ob- 
viously impossible to address specific lipids systems, but 
rather the more general properties of self-assembling sys- 
tems. 

The size of amphiphilic systems which can be simu- 
lated using simplified models is constantly growing by 
virtue of the availability of inexpensive and powerful com- 
modity PC hardware and due to the development of new 
simulation techniques such as Dissipative Particle Dy- 
namics (DPD) _23J. Simulations of model systems con- 
sisting of > 1000 lipids have been recently reported 
in the literature [23l 124 . |25| . The major restriction on 
the size of the systems in these simulation stems from 
the large number of atoms included in the simulation 
cell which is typically an order of magnitude larger than 
the number of amphiphilic molecules. The low ratio be- 
tween the number of lipids and the total number of atoms 
is due to two factors. The first one is the number of 
atoms comprising each lipid molecule which, as discussed 
above, can vary from 50 to 5 depending on the level of 
simplification employed in the simulations. The second 
factor is the number of water molecules in the simula- 
tion cell. In bilayers simulations the typical number ra- 
tio of water to lipid molecules is in the range of 10 to 
SOmilllllllliSEllil. A great fraction of the 
computing time is, thus, "wasted" on the simulations of 
the water even when the water molecules are represented 
by a single (unified) atom. Only very few models have 
so far been proposed in which the amphiphiles are sim- 
ulated without the presence of water. The major diffi- 
culty in establishing such "water-free" models is the need 
to mimic the hydrophobic effect that prevents the am- 
phiphilic molecules from leaving the aggregate into the 
solvent. Drouffe at al. and Noguchi at al. l27^ have 
used ad hoc multibody potentials to overcome this prob- 
lem. With the aid of these nonphysical potentials they 
have managed to observe the formation of fluid vesicles 
in their simulations. La Penna at al. have studied 
a water-free fiat bilayer model with Lennard-Jones (LJ) 
potentials that depend on the relative orientation of the 
lipids (and which are closely related to the Gay-Berne 
potentials used in liquid crystals simulations). With this 
model they have been able to simulate bilayer membranes 
in both the solid and the gel phases. Fluid membranes, 
however, were found unstable against lipids evaporation 
from the membrane plane. 

In this paper we present an exceptionally simple com- 
puter model of a fluid bilayer membrane. Our model has 
the following features: (a) It is a water-free model, i.e., 
we simulate the membrane without the presence of wa- 
ter, (b) The "lipids" forming the membranes consist of 
only three atoms, one representing the hydrophilic head- 
group and the other two the hydrophobic tail. These 
three atoms are "glued" to each other to form a rigid 
linear trimer (the lipid), and have no additional inter- 
actions between them, (c) The different lipids interact 
through finite range (truncated) LJ interactions between 



their three sites. The parameters of the LJ potentials 
are fixed and do not depend neither on the relative ori- 
entation of the lipids (as in Ref. .28]), nor on their local 
density (i.e., there are no multibody interactions in our 
model) . The above mentioned properties make our mem- 
brane model computationally very efficient (albeit a less 
"fiexible" one in comparison to other simplified models 
with more interaction sites per amphiphile). To inves- 
tigate the statistical mechanical properties of the mem- 
brane, we have performed a set of Monte Carlo (MC) 
simulations where for each MC run we have fixed the 
temperature, the number of lipids, and the projected area 
of the membrane. The projected area serves as the con- 
trol parameter in our simulations, and we have investi- 
gated the phase behavior of the membrane as a function 
of it. We found that upon increasing the projected area 
(i.e., reducing the area density of the lipids) the mem- 
brane undergoes a solid-fiuid phase transition. In the 
solid phase the lipids arc not mobile and they pack in a 
hexagonal order. In the fluid phase the lipids are free to 
diffuse in the membrane plane. We have measured the 
spectrum of thermal undulations of the fluid membranes 
from which we have extracted the surface tension and the 
bending modulus that characterize the elastic behavior of 
the membrane. At low area densities we found another 
transition from negative to positive surface tension, ac- 
companied by the formation of pores in the membranes. 
Such a behavior is indeed predicted by theoretical argu- 
ments jMEHIll 

The paper is organized as follows: In section ^ we 
present our computer model, and discuss the details of 
the simulations. In section ITTll we describe the physical 
properties of the systems as obtained by the simulations. 
The section is divided into three subsections dealing, re- 
spectively, with the phase diagram of the system, its spec- 
trum of thermal undulations and elastic properties, and 
the appearance of holes in membranes with large pro- 
jected area. We summarize and discuss the results in 
section Hvl 



II. DETAILS OF THE MODEL AND THE 
SIMULATIONS 

The lipids in our model system consist of three spher- 
ical atoms connected to form a linear trimer. The lipid 
molecules are rigid - they do not bend and the distance 
a between the center of the atoms is fixed (see fig. 
We set cr = 1 as our unit length scale throughout this 
paper. We shall label the three atoms forming each lipid 
as 1, 2, and 3. Atom 1 represents the hydrophilic head 
of the lipid, while atoms 2 and 3 represent its hydropho- 
bic tail. The different lipids interact with each other via 
spherically symmetric pair potentials between their con- 
stituting atoms. The pair potential Uij{r) depicts the 
interactions between atom i and atom j of two differ- 
ent molecules separated a distance r apart. The pair 
potentials U12 and Ui^ describe the interaction between 
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FIG. 1: A schematic picture of a lipid molecule in our model 
system - a rigid linear trimer consisting of three atoms whose 
centers are separated a distance a apart. The atom labeled 
1 (solid circle) represents the hydrophilic head of the lipid, 
while the atoms labeled 2 and 3 (open circles) represent the 
hydrophobic tail. 



hydrophobic and hydrophilic particles. They are given 
by the purely repulsive LJ potential 
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The pair potentials C/n, t/22, and t/33 describe the inter- 
actions between two similar atoms, both either hydropho- 
bic or hydrophilic. They are given by the attractive LJ 
potentials 
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where i = 1,2,3. Finally, the interaction between the 
hydrophobic particles 2 and 3 is also depicted by an at- 
tractive LJ potential, but of the form 



UM{r) = 46 
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(4) 



All pair potentials are truncated at the same cut-off 
Tc = 2.5(7, and the discontinuity at Vc is avoided by 
adding extra terms to the LJ potentials that ensure the 
vanishing of the potential, as well as of its first and sec- 
ond derivative, at r = r^. The final form of the pair 
potentials is thus given by 



Qj.2 



(r-rcf. (5) 



The different pair potentials are depicted in fig. 13 The 
values of the parameters (Jy (in units of a) and (in kT 
units, where T is the temperature and k the Boltzmann 
constant) used in the simulations are summarized in the 
caption on fig. [3 




FIG. 2: The different pair potentials Uij{r) used in our model: 
Uii (solid line), U22 (dashed line), U33 (dotted-dashed line), 
U12 (bold solid line), U13 (bold dashed line), and 1723 (bold 
dotted-dashed line). The distance r is in units of cr (see def- 
inition in text), while the potentials Uij are in kT units. 
The parameters aij and Sij [see Eqs.(0-(|^] are as follows: 

(711 = Llcr, (722 = 1.05(7, (733 = (7, (7l2 = L15(7, (7i3 = 1.4(7, 

a23 = 0.525(7, eii = 0.1875A;T, £22 = 1.75fcr, £33 = 1.875fcT, 
£12 = 1.1375fcT, £13 = 200fcT, and £23 = 375fcr. 



The pair potentials in our computer model have been 
designed to allow, on the one hand, the diffusion of 
molecules in the plane of the membrane but to restrict, 
on the other hand, their motion in the third direction. 
We have tested various models before we arrived to the 
one that we have used in the simulations. The origi- 
nal idea was to use dimers with one hydrophilic and one 
hydrophobic particles, and to describe the interactions 
between them by 6-12 LJ potentials [Eq.©] and a 12- 
power repulsive potential [Eq. (^] , depending on whether 
the atoms are of the same or different species. It turned 
out that the membranes depicted by such a model were 
unstable against the extraction of molecules from the 
membrane plane. To increase the membrane stability we 
added a third hydrophobic atom to the lipids. The pair 
interactions between this atom and the other two atoms 
are described by different forms of LJ potentials: For the 
interaction with the hydrophilic atom labeled 1 we use 
the more repulsive 18-power LJ potential U13 [Eq.®], 
while for the interaction with the hydrophobic atom la- 
beled 2 we use the 1-2 LJ potential U23 [Eq.®]. The 
former potential establishes a strong repulsion between 
the hydrophobic and the hydrophilic parts of the lipids, 
thus reducing significantly (eliminating on the time scale 
of the simulations) the escape probability of molecules 
(more on this point in the next paragraph). The latter 
has a very shallow minimum which allows a greater mo- 
bility of the lipids in the membrane plane (by making 
small the energy changes due to a relative motion of the 
lipids with respect to each other). We have gone through 
a rather lengthy "trial and error" process of fine tun- 
ing the parameters <7ij and eij which control the range 
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of pair repulsion and the depth of the attractive poten- 
tial wells. Their values have been set to (a) make favor- 
able the alignment of molecules next to each other at a 
distance slightly larger than cr, and (b) to make the at- 
traction between molecules sufficiently strong to support 
the stability of the membrane, but not too strong to the 
extent that would entirely prevent the diffusion of the 
lipids. 

It is not an easy task to form a fluid bilayer sheet 
in a model system that does not contain water. Mem- 
branes become fluid at low area densities and high tem- 
peratures, and under these conditions the lipids tend to 
escape quite easily from the membrane plane. It is the 
water that confines the lipids to the membrane. In the 
absence of water molecules this role has to be played by 
the hydrophobic heads which must form some sort of geo- 
metric or dynamic constraint for the extraction of lipids. 
In our model we establish such a constraint by making 
the excluded volume part of the pair potentials Uij non- 
additive, namely we make the size of a particle i "seen" 
by another particle i smaller than its size as seen by a par- 
ticle j of a different species. We can define the distance 
Uij at which the pair potential between them Uij = kT 
as a measure for the range of hard core repulsion between 
the two particles i and j. (It is unlikely to find a pair i 
and j separated by a smaller distance.) It is customary 
to regard the diameter of atom i and, with this 

interpretation, to expect for the additivity of the hard 
core diameters, i.e., to have a^- ~ {an -\- ajj)/2 for i ^ j. 
In our model we do not find this property (see fig. [21 . 
The pair potentials in our system describe the effective 
interactions between the different atoms and they include 
the effect of the water molecules which are not simulated 
explicitly. Therefore, there is no a priory reason why 
the effective diameters associated with different particles 
should be strictly additive. The increased range of hard 
core repulsion between the hydrophilic atom 1 and the 
hydrophobic atoms 2 and 3 is designed to compensate 
for the absence of water from the simulation cell. 

The simulations were performed with membranes con- 
sisting of = 1000 lipids (500 lipids in each layer) with 
periodic boundary conditions in the membrane (x, y) 
plane, and with no boundaries in the normal z-direction. 
Subsequent MC configurations were generated by two 
types of move attempts: translations of lipids and rota- 
tions around the mid (second) atom. The MC unit time 
is defined as the time (measured in number of MC config- 
urations) in which, on the average, we attempt to move 
and rotate each molecule once. The acceptance probabil- 



ity of both types of moves was approximately half. We 
performed a set of simulations of membranes with the 
same temperature T and number of lipids N , and with 
varying projected areas. For each value of the projected 
area we studied 8 different membranes starting at differ- 
ent initial configurations. The initial configurations were 
created by randomly placing 500 lipids in two layers with 
a vertical (along the normal z direction) separation a be- 
tween the atoms labeled 3 in the two layers, and with all 
the lipids oriented normal to the membrane plane. The 
initial configurations were then "thcrmalized" over a pe- 
riod of 5 • lO'^ MC time units, followed by a longer period 
of 6 • 10^ time units during which quantities of interest 
were evaluated. The duration of the MC runs is substan- 
tially larger than the relaxation time which we estimated 
in various ways: As a first approximation for the relax- 
ation time we used the time it took the potential energy 
of the membrane to saturate from its high initial value 
(resulting from overlap of particles in the random ini- 
tial configuration) to a final "typical" value. This time 
was of the order of 10"* MC time units. An indepen- 
dent estimate of the relaxation time was obtained from 
a study of the spectrum of thermal undulations of the 
membranes (see more details, later in the text). Inspec- 
tion of the autocorrelation function of the amplitude of 
the longest wavelength mode led to a similar estimate of 
10^ time units for the relaxation time. A more conser- 
vative estimate can be obtained from measurements of 
the self-diffusion constant of the lipids in the fluid phase 
(see, again, later in the text). The relaxation time can be 
associated with the time it takes a molecule to diffuse a 
distance equal to the pair potentials cut-off (2.5cr). The 
relaxation time obtained using this criterion was an order 
of magnitude larger (^ 10^ MC time units), still smaller 
than the equilibration time, and much smaller than the 
total length of the simulations. 



III. SIMULATION RESULTS 

A. Phase Diagram 

The projected area of the membranes in the simula- 
tions ranges from Ap = Ll = (26.875)^ to Ap = (30.625)^ 
with intervals of ALp = 0.625. For all area densities we 
measured the self-diffusion constant of the molecules rel- 
ative to the diffusion of the center of mass, defined by 



A '(fY 1 

D ^ hm -^-^ EE lim — ^ [{f,{t) - remit)) - (rl(0) - rcM(O))]' , (6) 

i—l 

where ri{t) denotes the position of the i-th lipid (defined by the position of its mid atom) at time t, while rcuit) 
denotes the position of the center of mass of lipids We have also measured the self-diffusion coefficient in the 
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membrane plane ^ defined by: 



1 ^ 

D^.y = l\m — J2 {[(^'^(t) - xcmit)) - (x^O) - a;cM(0))]' + [(y^it) - ycuit)) - (y.(0) - 2/cm(0))]'} , (7) 

where x and y denote Cartesian coordinates. (In all the simulations the membranes lied in the (x, y) plane, while 
fluctuating in the normal z-direction.) As the lipids can only diffuse within the plane of the membrane, we found no 
difference between D and Dx-v 
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FIG. 3: (a) Equilibrium configuration of a solid membrane 
with Ap = (26.875)^. Tlie atoms labeled 1 (the "hydrophilic" 
atoms) are shown as black spheres of diameter o, while the 
grey shaded spheres depict the (hydrophobic) atoms labeled 
2 and 3, and are of diameter a as well, (b) A top view of the 
plane of mid (labeled 2) atoms of the membrane upper layer. 

At low projected area (high area density) we found 
the membrane in a solid phase characterized by two fea- 
tures: (a) The diffusion constant of the lipids is vanish- 
ingly small. [The root mean square displacement V Ar'^ 
has barely changed during the course of the simulations, 
and it has never exceeded the typical distance between 
neighbor molecules a).] (b) The lipids arrange them- 
selves in a hexagonal order in the membrane plane. A 
typical equilibrium configuration of a membrane with 
Ap = (26.875)2 is shown in fig. (a). A top view of 
the plane of mid (labeled 2) atoms of the membrane up- 
per layer, revealing the hexagonal order of the lipids, is 
shown in fig. 01 (b). The lattice imperfections observed 
at fig. O (b) should be mainly attributed to the incom- 



mensurability of the 500 sites hexagonal lattice with the 
square simulation cell. 

At larger values of the projected area [Ap > (28.125)^] 
we found the membranes in a fluid phase. The main fea- 
ture that distinguishes fluid from solid membranes is the 
diffusion of the lipids. In flg.^we plot lipids mean square 
displacement Ar'^ [see definition in Eq.®] as a function 
of the simulation time t for fluid membranes with differ- 
ent area densities. The slope of the asymptotically linear 
curves is four times larger than the self-diffusion constant 
D. One can observe the growth of D with the increase 
of the projected area - a rather expected observation as 
the increase of the projected area means more room for 
the molecules to move. A typical equilibrium configura- 
tion of a fluid membrane with Ap = (28. 75)^ is depicted 
in fig. O (a) . Another characteristic feature of the fluid 
membranes is the loss of in-plane hexagonal order, as 
demonstrated in flg. El(b) [compare with flg. 131(b)]. 

The membranes with Ap = (30. 0)^ exhibited an inter- 
esting feature - they developed pores, as demonstrated 
in the configuration shown in fig. |^ These pores tended 
to appear irregularly in the membrane with a character- 
istic time scale r > 2 • 10^ for the formation of a pore, 
and a typical pore life time of a few thousand time units. 
Another interesting phenomenon which we observed for 
this value of Ap and did not observe at lower projected 
areas was the occurrence of "flip-flops" - the transition 
of lipids from one layer to the other. In flg. [71 we look 
at the same membrane depicted in fig. |S1 In this figure, 
however, we plot only the 500 lipids that were located 
in the upper layer in the initial configuration. About 30 
of them have managed to diffuse from the upper to the 
lower layer during the course of the simulations. A sim- 
ilar (although not necessarily identical) number of lipids 
have moved in the opposite direction. Trans-bilayer dif- 
fusion is an important process in real bilayer membranes 
[3^ . To allow for uniform bilayer growth, some of the 
lipids must be transfered from one leaflet to the other 
during the self-assembly process. When a flat bilayer is 
bent to form a spherical vesicle, the area of the inner layer 
becomes smaller than the area of the outer layer, and it 
is the transition of lipids from the former to the latter 
that balances their area densities. It has been suggested, 
based on experiments [s^ and computer simulations [3^ , 
that that the formation of pores and the flip-flop motion 
are closely interconnected. According to these studies, 
the pores provide a transverse diffusion conduit for the 
lipids, through which their hydrophilic headgroups cross 
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FIG. 4: The lipids mean square displacement /S.r'^ (measured 
in units) as a function of the time (measured in MC time 
units) for fluid membranes with (from bottom to top) A-p = 
(28.125)^ Ap = (28.75)^ Ap = (29.375)^ and Ap = (30.0)^ 



the hydrophobic region of the membrane. Our study sup- 
ports the conjecture about the possible relation between 
pores and flip-flops. In section lTlI Cl we discuss the origin 
of the formation of the pores. We show that it is associ- 
ated with the change in the sign of the membrane surface 
tension - from a negative to a positive value. Once the 
surface tension attains a sufficiently large positive value, 
the energy involved with the formation of the pores is 
compensated by the reduction in elastic energy. 

While the formation of the pores allowed the diffusion 
of molecules between the two layers, we did not observe 
[in membranes with Ap = (30.0)^] that pores also lead 
to the extraction of molecules from the membrane. It is, 
however, possible that the disassociation of pore- forming 
membranes occurs on time scales larger than the duration 
of our simulation. Fast disintegration of the membrane 
was observed when the projected area was increased to 
Ap — (30.625)^, which was, therefore, the largest pro- 
jected area set for the membranes in our study. 



B. Elasticity and Thermal Undulations 

On length scales larger than the membrane thickness, 
the bilayer can be modeled as a smooth continuous sheet. 
The thermal undulations of the bilayer can be studied 
with Helfrich Hamiltonian 9] relating the elastic energy 
to the shape of the membrane: 



dA 
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(8) 



The integration in the above equation is carried over the 
whole surface of the membrane. Three elastic moduli 
are involved with the Helfrich Hamiltonian: the surface 
tension 7, the bending modulus k, and the saddle-splay 
modulus kq. The quantities ci and C2 appearing in the 




FIG. 5: (a) Equilibrium configuration of a fluid membrane 
with Ap — (28.75)^. Black and grey atoms (of diameter a) 
depict hydrophilic (labeled 1) and hydrophobic (labeled 2 and 
3) atoms, respectively, (b) A top view of the plane of mid 
(labeled 2) atoms of the membrane upper layer. 



above equation are the local principle curvatures of the 
surface (see a rigorous definition in Ref.j^^l) which are 
surface invariants with respect to similarity transforma- 
tions (translations and rotations), while cq is the sponta- 
neous curvature of the surface. For flat bilayers cq = 0. 
It is customary to dispense with the use of the local cur- 
vatures in favor of two other (local) invariants: the mean 
curvature H = {ci + C2)/2, and the Gaussian curvature 
K = ciC2. If one only considers fluctuations which do 
not change the topology of the membrane, then the to- 
tal energy associated with the last term in Eq.(|SJl is a 
constantly. We, thus, arrive to the following more sim- 
plified form of Eq.®: 



/ (7 
Js 



(9) 



There are various ways to parameterize the surface. 
One of them is the Monge representation, where the sur- 
face is represented by a height function, z = h(x,y), 
above a reference x — y plane. For a nearly flat surface, 
i.e., when the derivatives of the height function with re- 
spect to X and y are small - h^, hy <^ 1, one obtains the 




FIG. 6: Equilibrium configuration of a fiuid membrane with 
Ap — (30.0)'^ having a pore on its upper right corner. 



following approximation for Eq. @ 
TC — dxdy 



^j{hl + hl) + ^. 



(10) 



Note that unlike Eq.®, the integral in Ea. HlOl) runs over 
the reference (x, y) surface rather than over the actual 
surface of the membrane. 

Equations (|5|)- ((Tn|l are expected to be valid only on 
length scales larger than the thickness of the membrane. 
The undulatory motion on smaller length scales (which 
we did not investigate in this study) is dominated by the 
so called "protrusion modes" js^ . In our simulations the 
profile of the bilayers was defined by mapping the system 
with linear size (of the projected area) L onto an 8 x 8 grid 
whose mesh size / = L/8 is indeed larger than the typical 
width of the membrane. The local height of the bilayer 
was then defined as the average of the local heights of 
the two layers. The latter were evaluated by the mean 
height of the lipids (whose positions were identified with 
the coordinates of their mid atoms) belonging to each 
layer, which were instantaneously located inside the local 
grid cell. The discretized form of Hamiltonian H1U|) is 



'H = 



(11) 



where summation goes over the discrete grid coordinates, 
and ao = P is the area of the grid cells. In Fourier 
coordinates we define 



and 



(12) 



(13) 



FIG. 7: Another view at the membrane depicted in fig. |^ 
Here we show only half of the lipids which were originally 
located on the upper leaflet of the bilayer. 



where the two-dimensional wavevector q has 8^ — 
64 discrete values satisfying {qx,qy = 2Trn/L, n — 
— 4, — 3, . . . , 2, 3}. In Fourier space the different modes 
decouple: 



(14) 



and, by invoking the equipartition theorem, we find that 
the mean square amplitude of the mode q (the "spectral 
intensity" ) 



kT 



(15) 



The instantaneous amplitudes of the different g-modes 
were evaluated using Ea. (|13ll once in every 100 MC time 
units, and were averaged over the course of the simu- 
lations. To extract the values of 7 and k we used the 
inverse from of Eq. (|15|) 
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n\q\*) 



kT 



(16) 



and plotted l/ao(|/if P) and a function of \q\'^. The re- 
sults of the spectral analysis of the undulations for fluid 
membranes with {Ap = 28.75)'^ (squares) and Ap = 
(29.375)^ (circles) are presented in fig.|Hl The error bars 
represent one standard deviation in the estimates of the 
averages. The curves depict the best fit of the numer- 
ical data to Ea. H16|) . obtained when 7 and k take the 
following values: 



7 



-1.4±0.2 



kT 



K = 54 ± 2 fcT, 



(17) 



for Ap = (28.75)2, and. 



kT 

7 = -0.6 ±0.2 — 
K = 42±2 kT, 



(18) 



for Ap = (29.375)2. We verified the vahdity of Eq.Jini) by 
attempting to fit our data to other polynomial functional 
forms, including a constants and a |g| ^ terms. The con- 
tributions of these terms to the fit were small, and did 
not result a significant change in our estimates of 7 and 
K, based on Ea. H16|l . 
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Iql^ [\I<J] 



The inverse of the spectral intensity for undulatory 

le 

(28.75)^ (squares) and Ap 



FIG. 

modes l/ao{\hq\^) as a function of the square wave number 



q for membranes with Ap 
(29.375)^ (circles). 



The above values of the bending modulus k are some- 
what larger than the values commonly reported in exper- 
iments in phospholipids: k ^ 10 — 20fcT [Toj. We should, 
in principle, correct our values of k by considering the re- 
duction in the effectively measured bending modulus due 
to long- wavelength thermal undulations. The correction 
term depends logarithmically on the size of the system 



elastic area energy more than they increase the elastic 
bending energy. For systems whose linear size is smaller 
than, but close to, Lc the amplitude of the small q modes 
become large [see Ea. (|15|l ]. and the approximation in 
Eas. Hll|l and H14() is no longer valid. It is then necessary 
to include higher order terms in the Hamiltonian H14|l . 
and to consider their influence on the spectrum of ther- 
mal undulations and on the stability of the membrane. 



C. Pore Formation 

The projected area Ap at which the surface tension 
vanishes is called the saturated {Schulman) area |40ll4ll|. 
One can evaluate the saturated area of our model mem- 
branes using a linear approximation for the relation be- 
tween the surface tension and the excess area 6Ap = 

A„-Ar,M 




(21) 



The coefficient Ka appearing in the above equations is 
the area compressibility modulus of the membrane. Us- 
ing the values of 7 obtained from the simulations for 
Ap = (28.75)2 ~ 826.6, and Ap = (29.375)2 ~ 862.9 
[Eqs.lEI) and (^HJ] in Eq.(|2U, we derive the following 
estimates: 



Ka = 19.6 ±6.6 



kT 



(22) 



K Attk V a 



(19) 



and, 



where a is some microscopic length. Setting a ^ I = L/8 
(the mesh size in our simulations), and using the values 
of K in Eas. l|17|) and (|18|l . we find that this correction 
amounts to about 1% of the value of k and, therefore, falls 
within the uncertainty in our estimates of the bending 
modulus. 

The fact that the surface tension 7 is negative has 
an interesting implication: It means that the size of the 
membrane cannot grow indefinitely, but an upper bound 
exists 



L, 



(20) 



above which there are small q modes with \q\ < qc = 
2tt/Lc that make the system unstable [see Ea. (|14f) ]. One 
can also understand the origin of this instability in "real" 
space, rather in q space: When 7 < the elastic energy 
of the membrane decreases by increasing its area, and it 
is the bending energy that stabilizes the system. Modes 
with larger wavelength (smaller \q\) require smaller cur- 
vatures to increase the area of the membrane and, thus, 
cost less bending energy. Membranes with linear size 
L > Lc have long wavelength modes that reduce the 



Ap = 890 ± 17, 



(23) 



for Ka and Ap. Equation H23|) suggests that our mem- 
brane with Ap — (30.0)2 _ QQQ ]-Qig}j^ ]-,(, found above 

the saturated area Ap and, therefore, may have a positive 
surface tension. The attempt to verify this conjecture by 
analyzing the spectrum of the membranes, as done for the 
fluid membranes with lower projected areas, is hampered 
by two technical difficulties. The first one is related to the 
fiip-flop motion of lipids between the two leaflets. The 
trans-bilayer diffusion makes it computationally compli- 
cated to determine which lipids are related to which layer 
of the membrane and, therefore, it becomes difficult to 
calculate the profile of the layers. The other difficulty re- 
sults from the holes which are created in the membrane. 
These pores may have an area larger than ag, the area 
of the grid cells. In such a case we find an empty cell 
with no lipids inside, and the height of the membrane at 
the corresponding grid point cannot the evaluated (un- 
less one interpolates this value using the height of the 
membrane at the adjacent grid points). We have taken 
advantage of the fact that the typical time for the ap- 
pearance of the pores and for the flip-flop motion which 
accompanies their formation, is larger than the relaxation 
time of the spectrum, and used short MC runs (during 
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which pores were not observed) to estimate the surface 
tension of the membrane. We found a positive surface 
tension with a magnitude of the order of 7 '--^ IfcT/cr^, 
which is roughly half an order of magnitude larger than 
the value anticipated by Eas. H21|l - H23|l . While the spec- 
tral analysis supports our conjecture that 7 is positive for 
Ap — (30.0)^, one should not attempt to use many inde- 
pendent short runs to achieve a more accurate estimate 
of 7. It is unclear how well equilibrated the membranes 
in these short MC runs are. Moreover, it is incorrect to 
base such an estimate on statistical averaging restricted 
to membranes without pores. The creation of the pores 
tends to reduce the surface tension since they make the 
effective area of the membrane smaller. 

Membranes with a positive surface tension can reduce 
their elastic energy by decreasing their area, and the for- 
mation of pores is obviously one of the mechanisms to 
achieve that. Other ways of reducing the membrane area 
which are not possible in our model is to decrease the pro- 
jected area or to increase the area density by adsorbing 
lipids from the solvent. For the case of a pore formation, 
one has to consider the line tension energy price involved 
with the creation of the hole. The simplest theoretical 
model discussing pore formation was suggested by Litster 
[Sfjj l . In this zero-temperature model, the contribution of 
a circular hole of radius i?porc to the free energy of a 
membrane with a positive surface tension 7 is given by 



F„. 



X2ttR^ 



(24) 



where A is the line tension of the hole. According to this 
model a pore with a radius larger than the critical value 
of A/7 is predicted to grow without bound. Such a ther- 
modynamically large circular hole can be created only if 
the critical energy barrier ttX^ /-f is accessible by thermal 
fluctuations. At a finite temperature it is necessary to 
take into account the entropy of the pores and the pic- 
ture becomes more complicated. Recent computer simu- 
lations have demonstrated that the typical shape of 
thermally induced pores is non-circular but rather of a 
self-avoiding ring or a branched polymer. The most strik- 
ing feature predicted by this study was, however, the fact 
that pores can appear at zero, and even at small negative 
surface tension. 

The major drawback of the above model is the fact 
that while it predicts the expansion of the pore without 
limit, the first term in Ea. (|24|l . assuming a linear relation 
between the reduction in elastic energy and the area of 
the pore, applies to small pores only. An improved model 
can be obtained by assuming other forms of the free en- 
ergy dependence on the pore area. We first consider a 
zero temperature model where the membrane does not 
fluctuate in the normal direction. The free energy of the 
membrane (which at zero temperature coincides with the 
potential energy) has a minimum at the saturated area 
which we shall now denote by [compare this notation 
with the one used in Ea. (|21|l ] to indicate that it is de- 
termined by energy consideration. The subscript p has 
been omitted since the projected area is also the total 



area of the membrane in this case. Close to A-^ we can 
use the quadratic approximation to describe the depen- 
dence of the free energy density / on the excess area 
SA^ = A~ A^ 



f 



F 



E 



fSA 



(25) 



where, as in the case of the saturated area, we use the 
superscript E in the notation of the area compressibility 
K^. If a pore of area Aporc is formed then the area of 
the membrane is reduced by Aporc and, consequently, the 
pore contribution to the free energy density is given by: 



fpOYC {A^ ApQYQ^ 



1 fSA^ -A 



pore 



V A^ 



2A0F 

A^ 



pore • 



(26) 



As in Eg. 1)24(1 . we consider a circular hole and, thus, its 
perimeter and area are related by F = vIttA. The equi- 
librium size of the pore A*^j.^ is found by solving the 
equation 9/poic/9Aporc = 0, and in addition by verify- 
ing that /pore ^;orc) < /pore 0) = 0. While in 
Litster's model a membrane with positive surface ten- 
sion can be only metastable against the formation of a 
pore, the model presented here yields a different scenario: 
Pores are thermodynamically unfavorable as long as the 
line tension satisfies 



A > A' = 



2 (M^)3/^ 
27^ Is 



(27) 



At this value a first order first transition occurs, and 
a pore of size Aporc — 2/3 SA^ is created. The pore 
grows gradually as A is decreased below this value. When 
A ^ 0, Apoic ~* SA^ , and the effective area of membrane 
attains the optimal (Schulman) value A^ . As in Litster's 
model, there exists a free energy barrier for the formation 
of the pore. At the transition (A = A') the height of the 
barrier is 



A4/3(M^)i/Vi^|'^'. 



(28) 



A theory for the entropic contribution to the free en- 
ergy of the pore has been recently presented by Sens and 
Safran ^A^. According to this theory, hole formation is 
one of the mechanisms to "redeem" the degrees of free- 
dom associated with the long wavelength modes in the 
fluctuation spectrum which are eliminated by the sur- 
face tension. To a first approximation, the effect of this 
entropic surface tension can be easily incorporated into 
Eqs.lEll) and Let us assume for a moment that 

= 0, and that the projected area of the membrane 
Ap is fixed. Because of the thermal fluctuations, the to- 
tal area of the membrane will be larger than Ap. As has 
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been explained in Ref. j4J|, there exists an optimal total 
area at which the membrane is tensionless 



A- 



All 



kT , / 











(29) 



where k is the bending modulus, and I is some molecular 
cut-off length. The superscript S denotes the fact that 
the optimal area discussed here is entropic in nature, and 
does not need to be equal to in Eq. H25(l . A membrane 
with A 7^ A^ will experience a surface tension. The free 
energy associated with this entropic surface tension can 
be calculated analytically. Here, however, we shall use 
the quadratic approximation in SA^ — A — A^ 



f 



F 



AT; - 



(30) 



which is valid only close to the minimum of the free en- 
ergy at A^ . The entropic area compressibility in Eg. 1301) 
is given by 0| 



(31) 



Combining the energetic H25|l and the entropic H30|) con- 
tributions to the free energy, we find another quadratic 
form for the total free energy of the membrane 



f{A) = — = -Ka (— 
Ap 2 y Ap 



(32) 



where the excess area 5 A = A — A is defined with respect 
to the minimum at 



^^K^Ap + KiA^ 
and the effective area compressibility is equal to 



A- 



(33) 



(34) 



The optimal area and the area compressibility appearing 
in the above two equations (and which include both en- 
ergetic and entropic contributions) should replace their 
purely energetic counterparts in Eq. H2t)|) for the pore free 
energy density and in Eas. (|27() and (|28|l for the criti- 
cal line tension and the free energy barrier. For typical 
values of phospholipids: k, — lOfcT ~ 5 x 10""'^'^ ergs, 
and Ap — (10 /im)^ = 10^^ cm^, we get upon substi- 
tution in Eq.JSU, ^ 5 x 10~^ ergs/cm^. This value 
of Kj^ is several orders of magnitude smaller than the 
area compressibility typically found in stretching exper- 
iments Ka ^ 10^ ergs/cm^ |l^ and, therefore, the en- 
tropic contribution to Ka and A can be neglected |45|. 
The effect of the thermal fluctuations can be felt only 
in small membranes (such as in this paper), as Kj^ in- 
creases by decreasing the projected area. Unfortunately, 



the membranes in our computer study are too small to 
allow the description of this effect by Eqs.lI2n|l-lI3IJ- The 
derivation of these equations is based on the assump- 
tion that the long wavelength behavior of the membrane 
is dominated by the surface tension. This requires that 
k{2ti/ yJ~A^Y ^ 7 - a criterion which is not satisfied in 
our case. The long wavelength fiuctuations in our mem- 
branes are mainly controlled by the curvature elasticity. 

Is the appearance of pores in the simulations in accord 
with the model described by Ea. H26(l ? In order to an- 
swer this question we need to evaluate the line tension 
A of the pore. The line tension A has the dimensions of 
energy per unit length. Its magnitude can be estimated 
by noting that the lipids on the rim of the pore have 
1-2 less neighbors compared to the other lipids. There- 
fore, the energy cost associated with each such lipid is of 
the order of the interaction energy between two adjacent 
molecules which is roughly kT. The length occupied by 
each lipid along the perimeter of the hole is of the order 
of (T, and so A ^ kT/a. This value of A should be smaller 
than the critical value A' given by Eq.imi). Using the 
values of Ka and A provided by Eas. iP^ and we 
arrive at the estimate A' ~ kT/a for the membrane with 



900 {5 A ~ 10). This means that A and A' are 



Ap = 30^ 

of the same order of magnitude and, thus, may obey the 
criterion given by Eq. H27|) for the thermodynamic stabil- 
ity of membranes with holes. The fact the pores in our 
simulations appear for only short time intervals, before 
they close up, may indicate that A is, in fact, slightly 
larger than A', and that the pores are only metastable. 
In addition to to the values of A and A', we also need 
to check the free energy barrier for the formation of the 
pores, as given by Ea. (|28|l . We find AF ^ kT, and so 
the opening of a pore can be nucleated by thermal fluc- 
tuations. 



IV. SUMMARY AND DISCUSSION 

We have introduced a new simple computer model for 
bilayer membranes whose main feature is the fact that the 
system is simulated in vacuum rather than in aqueous en- 
vironment. The elimination of the solvent from the simu- 
lations greatly improves computational efficiency. Devis- 
ing a "water-free" model is a great challenge as the wa- 
ter molecules, via the electrostatic interactions between 
them and the lipids, play a central role in the aggregation 
and the stabilization of the membrane through the result- 
ing hydrophobic effects. The self-assembly of the system 
has not been investigated in this paper. (The reader is 
referred to the simulations presented in Refs. [2^ and 
in which this issue has been addressed). We did, how- 
ever, demonstrated that bilayers, once they are formed, 
can be stable without the surrounding solvent. One only 
needs to modify the interactions between the lipids, and 
use effective potentials that compensate for the absence 
of water by producing a barrier against the disintegra- 
tion of the membrane. In this model we have been able 
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to stabilize the membrane using pair-wise short range 
interactions only - another feature that reduces the com- 
putational effort involved with the simulations. To the 
best of our knowledge, this is the first water-free com- 
puter model in which fluid membranes are being observed 
without the need of multibody interactions. 

We have found that our simple model reproduces many 
known features of bilayer membranes, such as the tran- 
sition from a high density solid phase to a low density 
fluid phase. We have inspected the spectrum of thermal 
fluctuations of the fluid membranes, and found it to be 
well described by the Helfrich Hamiltonian. From the 
analysis of the spectral intensity of the different modes, 
we have extracted the surface tension and bending mod- 
ulus of the system. Based on our numerical results for 
the surface tension, we have attempted to determine the 
optimal area of the membrane at which the surface ten- 
sion vanishes. Indeed, for areas larger than the saturated 
area, we have found evidences that the surface tension be- 
comes positive. Fluid membranes with positive surface 
tension can develop pores, and the creation of pores al- 
lows the diffusion of lipids from one layer to the other 
(flip-flops) . The opening of holes in our membranes is in 
agreement with a simple model that takes into account 



the contributions to the area compressibility of both the 
inter-particle forces and the thermal fluctuations. 

In order to make a closer contact with biological sys- 
tems, it is necessary to extend the model presented here 
to include the other elements found in biomemebranes 
such as the membranes proteins and the cytoskeleton. 
It would be interesting to see whether these additional 
components can also be modeled in a coarse-grain man- 
ner that would minimize both the computational and the 
conceptual complexity. Such a model may shed light on 
an abundance of challenging problems like the effect of 
the cytoskeleton on the elastic properties of the bilayer, 
or the role played by the membrane proteins in transport 
processes across the membrane. 
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